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ABSTRACT 

New subgrid-scale models for the large-eddy simulation of compressible turbulent flows are 
developed and tested based on the Favre-filtered equations of motion for an ideal gas. A 
compressible generalization of the linear combination of the Smagorinsky model and scale- 
similarity model, in terms of Favre-filtered fields, is obtained for the subgrid-scale stress ten- 
sor. An analogous thermal linear combination model is also developed for the subgrid-scale 
heat flux vector. The two dimensionless constants associated with these subgrid-scale models 
are obtained by correlating with the results of direct numerical simulations of compressible 
isotropic turbulence performed on a 96 3 grid using Fourier collocation methods. Extensive 
comparisons between the direct and modeled subgrid-scale fields are provided in order to val- 
idate the models. A large-eddy simulation of the decay of compressible isotropic turbulence 
- conducted on a coarse 32 3 grid - is shown to yield results that are in excellent agreement 
with the fine grid direct simulation. Future applications of these compressible subgrid-scale 
models to the large-eddy simulation of more complex supersonic flows are discussed briefly. 
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1. Introduction 


The direct numerical simulation of turbulent flows at the high Reynolds numbers encountered 
in problems of technological importance is all but impossible as a result of the wide range 
of scales that are present. Consequently, the solutions to such problems must invariably 
be based on some form of turbulence modeling. Traditional turbulence models based on 
Reynolds averages have had only limited success since the large scales of the turbulence 
- which contain most of the energy - are highly dependent on the geometry of the flow 
being considered. Experience has indicated that such models usually break down when a 
variety of turbulent flows are considered (Lumley 1983). The small scales are more universal 
in character, and serve mainly as a source for dissipation. Hence, it can be argued that a 
better understanding of turbulent flows could be achieved if just the small scales are modeled 
while the large scales are calculated (Deardorff 1970). This is the fundamental idea behind 
large-eddy simulations. 

During the past decade, considerable progress has been made in the large-eddy simulation 
of incompressible turbulent flows. This effort has shed new light on the physics of turbulence. 
The earliest work relied heavily on the use of the Reynolds averaging assumption to elim- 
inate the Leonard and cross stresses while the Reynolds stresses were computed using the 
Smagorinsky model (Deardorff 1970, Leonard 1974, Reynolds 1976). More recent large-eddy 
simulations have been based on the direct calculation of the Leonard stresses with models 
provided for the cross and Reynolds subgrid-scale stresses in order to enhance the numerical 
accuracy (see Biringen and Reynolds 1981, Bardina Ferziger and Reynolds 1983). However, 
among these newer models, only the Bardina, Ferziger and Reynolds (1983) model, with a 
Bardina constant of 1.0, satisfies the important physical constraint of Galilean invariance 
(Speziale 1985). The underlying physical concepts, fundamental numerical algorithms, and 
comprehensive historical data behind the recent field of large-eddy simulation have been 
presented in articles by Schumann (1975), Voke and Collins (1983) and Rogallo and Moin 
(1984). More recently, work on the subgrid-scale modeling of transition to turbulence of ini- 
tially laminar incompressible flows has begun (Piomelli, Zang, Speziale and Hussaini 1990). 
Several large-eddy simulations have been performed and initial results are promising. 

Despite the intensive research effort that has been devoted to the large-eddy simulation 
of incompressible flows as outlined above, it appears that no large-eddy simulation of a 
compressible turbulent flow has yet been attempted. Of course, such work could have im- 
portant technological applications in the analysis of turbulent supersonic flows, where shock 
waves are generated, and in turbulent flows within combustion chambers. The prerequisite 
for carrying out such computations is the development of suitable subgrid-scale models for 
compressible turbulent flows. With the exception of the recent work of Yoshizawa (1986) 
and Speziale et al. (1988), few, if any, studies along these lines appear to have been pub- 
lished. The subgrid-scale models of Yoshizawa are only suitable for slightly compressible 
turbulent flows since they made use of an asymptotic expansion about an incompressible 
state. Recently however, Dahlburg, Zang and Dahlburg (1990) have performed an extensive 
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parameter study using the model developed by Speziale et al. (1988). It was shown by 
Speziale et al. (1988) that, for the purposes of accuracy, the Leonard and cross stresses 
must be accounted for. Furthermore, the modeling of the isotropic part of the Reynolds 
subgrid-scale stress tensor was shown to be questionable - an issue that was left for future 
research. 

In this paper, complete subgrid-scale models are developed for the closure of the Favre- 
filtered Navier-Stokes and energy equations. The compressible subgrid-scale stress model 
that is obtained in Section 2 reduces to the linear combination model of Bardina, Ferziger 
and Reynolds (1983) in the incompressible limit. Likewise, the subgrid-scale heat flux model 
that is obtained herein consists of an analogous linear combination of scale similarity and 
gradient transport terms. The dimensionless constant which appears in the subgrid-scale 
stress model is arrived at through correlation analysis of data generated from direct numerical 
simulations of compressible isotropic turbulence. A more detailed comparison of computed 
and modeled subgrid-scale fields is presented along with the results of a large-eddy simulation 
of compressible isotropic turbulence. 


2. Subgrid-Scale Models for Compressible Turbulence 

The compressible turbulent flow of an ideal gas is considered. Such flows are governed by 
the continuity, momentum and energy equations which - neglecting body forces - are given 
by (cf. Batchelor 1967) 

d P . d{pv k ) 

dt dx k 


d(pvk) d(pvkvi ) 


dt 

d(ph ) d(phv k ) 


dxi 


dp dajci 
dx k dxi 


( 1 ) 

( 2 ) 


dp dp 

= dt + Vk d V k 


d , dT , 

+ ~ — (« « — ) + $ 
dx k dx k 


( 3 ) 


dt ' dx k 

respectively, where p is the mass density, v is the velocity vector, p is the thermodynamic 
pressure, p is the dynamic viscosity, h is the enthalpy, T is the absolute temperature, and k 
is the thermal conductivity. The viscous stress a k i and the viscous dissipation $ are defined 

by 
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respectively. Herein, the Einstein summation convention applies to repeated indices. Equa- 
tions (l)-(3) must be supplemented with the equations of state 


P = pRT, 


C P T 


( 6 ) 
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for an ideal gas where R is the ideal gas constant and C p is the specific heat at constant 
pressure. Likewise, the dependence of the viscosity and thermal conductivity on the temper- 
ature must be provided (i.e., relationships of the form /i = fi(T ) and k = k (T) are needed 
and these depend on the gas under consideration). 

Any flow variable can be filtered in the following manner. 

T{x) = [ G(x — z,A)^F(z)d 3 z CO 

where G is a filter function, A is the computational mesh size, and D is the domain of the 
fluid. The filter function G is typically taken to be an infinitely differentiable function of 
bounded support in a bounded domain, or a Gaussian distribution in a periodic domain . It 
is normalized by requiring that 

f G(x — z, A)d 3 z = 1. W 

J D 

It follows that in the limit as the computational mesh size goes to zero, (7) becomes a Dirac 
delta sequence, i.e. 

lim / G(x - z, A)F(z)d 3 z = [ 6(x- z)F(z)d 3 z = F(x) (9) 

A— *-0 Jo J D 

where S(x-z) is the Dirac delta function (Arfken 1970). The filter function has the property 
that the amplitude of the high-frequency spatial Fourier components of any flow variable T 
are substantially reduced. Consequently, T represents the large-scale part of T. At this 
point, it should be mentioned that as a result of the defining properties of G, it follows that 

d]F_ dT / 10) 

dt dt ’ dxk dxk 

Piomelli, Ferziger and Moin (1987) discuss the relationship between the form of the filter 
function and that of the subgrid-scale turbulence model. 

The turbulent fields are decomposed as follows, based on Favre filtering: 

+ (11) 


where the Favre filter 



( 12 ) 


is defined in an analogous manner to the Favre time average which has been of use in the 
more traditional studies of compressible turbulent flows (Hinze 1975). However, contrary to 
the more traditional Favre time averaging, 


t A Gaussian filter is adopted for the calculations in this study 


( 13 ) 
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in general, and hence 

r ± o. (i4) 

The direct filtering of the continuity equation (1) yields 

where we have used (10) and (12). Likewise, a direct filtering of the momentum equation 
yields 

d(pvk) d(pv k vi) dp da kl dr kl 

dt dxi dx k dxi dxi 


where 


p = pRT 


Tu = ~p(vkVi - v k vi + v' k vi + v\v k + v' k vl) (18) 

is the subgrid-scale stress tensor. The subgrid-scale stress tensor can be decomposed as 
follows, 

t = L -f* C -f- R (^-®) 

where 

Lu = -p{v k vi - v k vi ) (20) 

Ci u = —p(v k vi + v\v k ) (21) 

Rki = -pv' k v[ ( 22 ) 

are respectively, the subgrid-scale Leonard, cross, and Reynolds stresses based on Favre 
filtering. From (20), it is clear that the Leonard stress can be calculated directly and does 
not need to be modeled. The cross stress is modeled with the scale similarity model 

C k i = -p(v k vi - v k vi ) (23) 

(with a coefficient of unity to ensure Galilean invariance of the overall model). This model 
is analogous to its incompressible counterpart, which has been reasonably successful in the 
large-eddy simulation of incompressible turbulent flows (Bardina, Ferziger and Reynolds 
1983, Speziale 1985). The subgrid-scale Reynolds stress tensor is separated into deviatoric 
and isotropic parts, respectively, as follows: 

R = r>R + /R (24) 


where 


D R k[ = -p{v' k v\ - -v-v-Ski), 

iRki = —pv'iV'iSki. 
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Here, the deviatoric part of the subgrid-scale Reynolds stress tensor, dR, is modeled using 
the compressible generalization of the Smagorinsky model that is given by 


pRki = 2C R pA 2 II 1 § /2 (S kl - jrS mm 8 k t) 


1 - 


where 


c _ l t d ' 3 * , d ^\ 

Sh 2^dxt + dxj 


(27) 

(28) 

n$ = S mn S mn (29) 

(i.e., S is the Favre filtered rate of strain tensor while II§ is its second invariant) and Cr is 
the compressible Smagorinsky constant. Yoshizawa - by means of a two-scale DIA method - 
derived a model for the isotropic part of the subgrid-scale Reynolds stress tensor, j R, given 

by 2 

iRu = -3 CjpA 2 II s 8 kl (30) 

where Cj is a dimensionless constant. Equation (30) can, for the most part, be obtained 
by making a turbulence production equals dissipation equilibrium hypothesis (Yoshizawa 
1986). However, this model was shown by Speziale et al. (1988) to correlate very poorly 
with the results of direct numerical simulations of compressible isotropic turbulence. Since 
jR is extremely small compared to the thermodynamic pressure, we propose to neglect it - 
an assumption that will be justified later. Hence, the overall subgrid-scale stress model we 
propose takes the form 


r k i = ~p{v k vi - v k vi) + 2C R pA 2 II 1 s /2 (S k i - | S mm , 

In the incompressible limit, Eq, (31) reduces to the linear combination model 

- T W /p = m~%vi - 2 C r A 2 IIj 2 S k i 


1 - 


M- 


(31) 
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of Bardina, Ferziger, and Reynolds (1983) where the Bardina constant is one in order to sat- 
isfy Galilean invariance (Speziale 1985). This reduction process is a consequence of Eq. (31) 
and 

V = V, Smm = S'mm = 0 (33) 

when p becomes constant. 

A direct filtering of the energy equation yields the filtered form 
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is the subgrid-scale heat flux. The subgrid-scale heat flux can be decomposed in the same 
fashion as the subgrid-scale stresses. This leads to 


where 


Q = QW + Q( c ) + Q(*> 

(36) 

Q (l) = C p p{W?-vf) 

(37) 

Q( c > = c p p{v' k f + iZr>) 

(38) 

Q (H) = C p pv^T' 

(39) 


are the Leonard, cross, and Reynolds heat fluxes. Analogous to the modeling of the cross 
stress, the cross heat flux is modeled using the scale similarity format 


Qk C) = C p p(v k T - h f). 


(40) 


The Reynolds heat flux is modeled with the usual gradient transport format as follows (cf. 
Eidson 1985): 

<«> 

where Pr t is the turbulent Prandtl number. Of course, the Leonard heat flux can be 
calculated directly. Hence, the overall model for the subgrid-scale heat flux we propose is as 
follows: 


Qk — C pP 


(v k T - v k T) - 
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Pr T s dx k 
and is obtained by combining equations (37), (40), and (41). 


(42) 


At this point, some comments need to be made concerning the viscous terms on the 
right-hand side of (16) and the pressure gradient- velocity and viscous dissipation terms 
which appear on the right-hand side of (34). The pressure gradient-velocity correlation can 
be written in the alternative form 
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,fK 

dx k 

model and not much success has been achieved in dealing with it in the context of Reynolds 
stress models. However, within the framework of subgrid-scale modeling, this term and its 
corresponding cross correlation have physical interpretations. They represent the contribu- 
tion of the dilatation of the small scales to the internal energy variation of the fluid - an 
effect which is expected to be small. Hence, for this initial study, we neglect such terms. 
Furthermore, since the mean temperature is constant and the temperature fluctuations are 
small (< 10%), the viscosity and thermal conductivity are held constant. For similar reasons, 
we also neglect the small scale component of the viscous dissipation. 

The turbulence model proposed herein is thus complete once values for the constants 
Cr and Pry are obtained. This will be accomplished using the results of direct numerical 
simulations of compressible isotropic turbulence. 


is not yet closed. The temperature dilatation correlation (T 


) is extremely difficult to 


3. Numerical Method 


Our direct simulations of compressible turbulence are based on a non-dimensional form of 
Eqs. (l)-(3), with the time derivative in the energy equation written solely in terms of the 
pressure. In order to alleviate the severe stability limit imposed at very low Mach numbers 
by the acoustic waves, a splitting method is adopted. The first step integrates the equations 


a JL =0 

dt 

(45) 

d(pv Jfc) d(pv k vi ) _ da kl 

dt dxi dxi } 

(46) 

dp \ v k dp 1 ~<v dVk 1 9271 1 r T rn 

dt dxk dx k 0 dx k RePrM ^ dx^dx^ ’ 

(47) 

while the second step integrates 

dp d(pv k ) 
dt dx k 

(48) 

d(pv k ) dp 

dt dx k ’ 

(49) 

dp , _2 d {P v k) _ n 
dt 0 dx k 

(50) 


The constants c 0 and M ^ are the current root mean square (rms) value of the sound speed (c) 
and the reference Mach number, while 7 = C p /C vy where C v is the specific heat at constant 
volume. These equations are non-dimensionalized in terms of a length scale (L 0 ), a velocity 
scale (U 0 ), a pressure scale (P 0 ), a reference viscosity \i and a reference thermal conductivity 
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k. The Reynolds number is given by Re = p 0 U 0 L 0 /fi, the Prandtl number by Pr = C p p/k, 
and the reference Mach number by M , «, = U 0 I{^RT 0 ) x/2 . For all calculations presented in 
this study, 7 = 1.4, and Pr = 0.7. Initially, the density p 0 is uniform and equal to one. The 
computational domain is a cube, normalized to [0,27 t] 3 . Periodic boundary conditions are 
imposed in all three directions. 


The spatial derivatives in these equations are approximated by a Fourier collocation 
method (see, for example, Hussaini and Zang 1987). In each coordinate direction, N grid 
points are used: x kj = 2 tt j/N, for j = 0, 1, ..., N — 1. The derivative of a function ^(x) with 
respect to x k is approximated by the analytic derivative of the trigonometric interpolant of 
^(x) in the direction x k . Most simulations of incompressible, homogeneous turbulence have 
used a Fourier Galerkin method. The compressible equations, however, contain cubic rather 
than quadratic nonlinearities and true Galerkin methods are more expensive (compared with 
collocation methods) than they are for incompressible flow. The essential difference between 
collocation and Galerkin methods is that the former are subject to both truncation and 
aliasing errors, whereas the latter have only truncation errors. As discussed extensively by 
Canuto, et al (1987), the aliasing terms are not significant for a well-resolved flow. However, 
care is needed to pose a collocation method in a form which ensures numerical stability. For 
this reason, the second term in (46) is actually used in the equivalent form 


1 

2 


d{pv k vi ) dv k 

~“&T“ + pv ‘d7, 


+ Vk 


djpvi) 

dxi 


(51) 


As noted by Feiereisen, Reynolds and Ferziger (1981), when this form is employed together 
with a symmetric differencing method in space (for example Fourier collocation), then in 
addition to mass, and momentum, energy is conserved for the ideal compressible equations 
(zero viscosity and thermal conductivity) in the absence of time differencing (and splitting) 
errors. 


The second fractional step of the splitting, given by (48)-(50), contains most of the effects 
of the acoustic waves. This splitting is employed at each stage of a third-order Runge- 
Kutta method. In the simulations reported here, the second fractional step is integrated 
analytically. In Fourier space, (48)-(50) become 


dp 

+ tkimi = 0 


(52) 


drhi 

~dt 


+ ikip = 0 


^ + ic 2 0 kirhi = 0 


(53) 

(54) 


where mi = pvi and Fourier transformed quantities (which depend upon the wavenumber k) 
are denoted by a circumflex. The exact solution of these equations is 



cos(cokAt a ) + B sin(co&At,) — A) 


(55) 
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( 56 ) 


rk\ 2 ^ — m\^ — '~r[A sin(coA:At J ) — B cos(c 0 kAt 8 ) + B ] 
c 0 k 

p( 2 ) = A cos(c 0 kAt B ) + B sin(c 0 Ai,) (57) 

where k = |k| is the magnitude of the Fourier wavenumber, A = p^\ B = The 

superscript 1 denotes the result of the first fractional step of the splitting and the superscript 
2 the results of the second fractional step. The effective time-step of the Runge-Kutta stage 
is denoted by At,. The advantage of this splitting is that the principal terms responsible for 
the acoustic waves have been isolated. Since they are treated semi-implicitly, one expects 
the time-step limitation to depend upon v + \c — Co\ rather than v + c. (Although there is 
also a viscous stability limit for the first fractional step, it is well below the advection limit 
in the cases of interest.) This is clearly a substantial advantage at low Mach numbers. Since 
the second fractional step is integrated analytically, it does not contribute to any time step 
limitations. If one is truly interested in all the details arising from the sound waves, or if 
there is a substantial coupling between the sound waves and the rest of the flow, then the 
time-step must be small enough to resolve the temporal evolution of these waves. But, if 
only the larger-scale sound waves are of interest, then this splitting method is useful. 

During the acoustic fractional step, an isotropic truncation is performed: for each variable 
(p, pv and p), all Fourier coefficients for which 

kik t > (AT/2) 2 (58) 

are set to zero. This reduces the numerical anisotropy produced by a cubic truncation. 
Moreover, it reduces the aliasing interactions in the collocation method (Canuto et al 1987, 
Chapters 3 and 7). 

The compressible code can also be executed in a purely explicit mode. In this case no 
splitting is performed; Eqs. (45)-(50) are simply combined in the appropriate manner and 
integrated directly. 


The expected stability limit of this three-dimensional Fourier collocation method for the 
compressible Navier-Stokes equations has the form 


A ‘ < ° [■£? g w '] (f) 

(59) 

where for the semi-implicit version, 


Ui = |u,| + |c - c 0 | 

(60) 

while 


Ui = |«i| + |c| 

(61) 


when the time advancement is fully explicit. For the third-order Runge-Kutta method 
employed here, we use a = 0.5. 


A number of simulations have also been conducted of strictly incompressible flow. These 
were performed with a separate code which also used a Fourier collocation method, but for 
the simpler, incompressible Navier-Stokes equations. 
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4. Comparison with Incompressible Results 


The initial conditions for the numerical simulations were designed to reproduce the experi- 
mental data of Comte- Bellot and Corrsin (1971) on isotropic turbulence, hereafter referred to 
as CBC. These experiments were also the basis of direct simulations used by Clark, Ferziger 
and Reynolds (1979), Bardina, Ferziger and Reynolds (1983), and McMillan and Ferziger 
(1979) in their analyses of incompressible LES models. Initial conditions are chosen to match 
CBC measurements at a non-dimensional time of 240 (cf. table 4 in Comte- Bellot and Corrsin 
1971). The computational domain is a cube of side 20/27T cm. The CBC parameters are 
associated with measurements taken behind a grid with a mesh spacing of one inch, and a 
mean fluid velocity of 393.7 in/ sec. The initial time in the direct simulation corresponds to 
t = 0.00254 sec in the CBC experiment. The reference length {Lq), velocity (Uo) and pres- 
sure (P 0 ) are respectively 20/ 27T cm, 1 cm/ sec and 1 gr j am sec 2 . After generating a random, 
divergence-free velocity profile, the kinetic energy (in Fourier space) is scaled to match the 
measured CBC energy spectrum (see Appendix). Finally, the velocities are scaled so that 
the initial rms velocity agrees with the measured values. This adjustment is typically less 
than 1%, which provides one measure of the uncertainty in the fit to the experimental data. 
With the chosen non-dimensionalization, the Reynolds number Re = UqLq/v is 22.74 based 
on a kinematic viscosity v — 0.14 cm 2 /sec. Table 1 summarizes the parameters measured 



CBC 

64 3 

~96 t ~ 

1W 

Vrmi 

6.75 

6.75 

6.75 

6.75 

\tr(Sk) 

- 

0.0 

0.0 

0.0 

E 


68.3 

68.3 

68.3 

e 

462 

375 

432 

447 


0.26 

0.28 

0.27 

0.25 

A i2 

- 

0.20 

0.19 

0.27 

■^13 

- 

0.20 

0.19 

0.26 

R A 

38.1 

43.6 

41.3 

37.8 


Table 1: Initial conditions based on CBC experiment and Clark et al. 
(1979) calculation. Mach number is zero. 


by CBC at t=240. The Taylor microscale length A y is defined by 



( 62 ) 
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and the dissipation e by 


e = 2 n J SijS{jd 3 x, 

(63) 

where Sij is the rate of strain tensor 


1 (dvi dvA 
ij ~ 2 ^ + dxj ' 

(64) 

In Eq. (62), (•) denotes a spatial average. Its exact definition is given in 
Taylor microscale Reynolds number is 

the Appendix. The 

R> = V ' Xu - 

(66) 


V 


The velocity derivative skewness and flatness tensors Sk and FI are the third and fourth 
moments of the velocity gradient and are defined by 



In tables 1 and 2, only the trace of the skewness and flatness tensors are shown. The 
remaining columns list the parameters obtained from the initial conditions of the numerical 
simulations on 64 3 , 96 3 and 128 3 grids. There is a 20% discrepancy between the dissipation 
obtained by CBC and the dissipation computed on the coarsest grid which suggests that a 64 3 
grid has marginal resolution, at best. A 12% difference between the value of R\ calculated 
on the 64 3 grid and that obtained by CBC confirms the need for grids finer than 64 3 . On a 
96 3 grid, both the dissipation and R\ are in much closer agreement with CBC. Discrepancies 
between our results and CBC for e and Rx are respectively 6.5% and 7.5% on a 96 3 grid. On 
the finest grids on which the direct simulations were performed, the computed values of e 
and R\ , respectively, have relative errors of 3.5% and less than 1% when compared to CBC. 

The numerical simulations were run from t — 240 until t = 375 (in CBC units), which 
corresponds to a non-dimensional time interval of 0.1145 (in our units). Table 2 furnishes a 
comparison of the experimentally measured parameters with those from the numerical sim- 
ulation at the final time. On the coarsest grid, the total dissipation rate that was calculated 
is still slightly below the value measured by CBC. A 96 3 grid generates values of e consistent 
with CBC. 
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CBC 

64 3 

96 3 

128 3 

^rmj 

5.03 

5.18 

5.19 

5.21 

IM Sk) 

- 

-0.42 

-0.51 

-0.52 

E 

38.6 

40.3 

40.4 

40.7 

e 

154.4 

151.3 

154.6 

156.8 


0.34 

0.33 

0.34 

0.33 

^12 

- 

0.24 

0.24 

0.23 

^13 

- 

0.24 

0.24 

0.23 

Rx 

36.6 

37.8 

40.4 

37.2 


Table 2: Final conditions (t=0.1145) based on CBC experiment and 
Clark et al. (1979) calculation. Mach number is zero. 


At t=0.1145, the diagonal components of Sk are —0.5 which agree well with the numerical 
results of Kerr (1985). Kerr studied isotropic, turbulent flow, but prevented the decay of 
energy by using an exterior energy source at the large length scales. 

As noted in the previous section, we have chosen not to de-alias the advection terms. 
In reaching this decision we drew upon the extensive evidence that has accumulated on 
aliasing effects in the last dozen years (Canuto, et al 1987., Chapters 3, 4 and 7) and upon 
tests conducted with the incompressible isotropic turbulence code. In this code, de-aliasing is 
accomplished by applying the 2/3-rule (Canuto, et al 1987, Chapters 3 and 7) in an isotropic 
fashion; e.g., the de-aliased results for a 64 3 grid are obtained by running the incompressible 
code on a 96 3 grid and applying the truncation given by (58) with N/3 in place of N/2 on 
the right hand side. The results are summarized in Fig. 1. Here we present the energy 
spectra E(k ) (defined in the Appendix) for 64 3 , 96 3 , and 128 3 grids at t = 0.0586 for both 
aliased and de-aliased calculations. Some adverse effects of aliasing are apparent on the 64 3 
grid, but they are only in the tail of the spectra, and they are already insignificant on a 96 3 
grid. For the reasons outlined here, a 96 3 grid was chosen as the standard discretization for 
the incompressible and for the compressible simulations. 


5. Compressible Turbulence Results 

5.1. Direct Simulations 

Recent work on the direct numerical simulation of homogeneous compressible turbulence 
has indicated the crucial role played by the initial conditions. Passot and Pouquet (1986) 
conducted direct simulations of two-dimensional, compressible isotropic turbulence and con- 
cluded that when the initial rms density fluctuations are small, the turbulence statistics 
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remain quasi-incompressible for turbulent Mach numbers Mt less than 0.3. 

M t = (^-J . (68) 

They also demonstrated (through the use of direct numerical simulations) that eddy shock- 
lets result for sufficiently high initial rms density fluctuations and/or turbulent Mach num- 
bers. A more systematic analysis and categorization of the effect of the initial conditions on 
compressible isotropic turbulence was achieved recently by Erlebacher et al. (1990). 

They concluded that for 0 < Mt < 0.3, p T m§ must initially be 0{Mt) for the resulting 
turbulence statistics to become strongly compressible with an 0(1) ratio of compressible 
to incompressible turbulent kinetic energy. (For the range of Mt considered herein, no 
eddy-shocklets occur.) On the other hand, if initially p T mti ^rmi ^ Mt, then the resulting 
turbulence statistics remain quasi-incompressible. 

We first present the results of direct numerical simulations of compressible isotropic 
turbulence corresponding to the initial conditions of the CBC experiment but with a variety 
of non-zero mean Mach numbers. Since for these simulations, the initial conditions are 
p Tma = 0, T Tma < 1, only weakly compressible turbulence statistics are expected according 
to the theoretical results of Erlebacher et al. (1990). Unless specified otherwise, a subscript 
rms for any variable T refers to the quantity [{T - (•7 r )) 2 ) 1/3 /(«7 r )- 

The initial pressure distribution over the entire field is specified. The fluctuating pressure, 
p f is determined from the velocity distribution by enforcing a zero initial time derivative for 
V • v. A Poisson equation for p f is obtained from the divergence of the momentum equation 
after setting the time variation of V • v to zero (Feiereisen et al. 1981). The mean pressure, 
p TO , is then determined so that a prescribed initial mean average Mach number, M 0 , defined 
to be the ratio of rms fluid velocity and rms speed of sound, is achieved. An analytic 
expression for p m is given by 

_ Mp l v 2 d 3 x + / 7 P//>~ x d 3 x / 6 g\ 

Pm ~ 7 fp-Wx 

The initial average Mach number is specified at the outset of the direct numerical simulations 
(DNS) as an initial condition. Density is initially set to unity, while the temperature, if 
required, is derived from the equation of state. Direct numerical simulations are performed 
for M 0 = 0.0, 0.1, 0.4 and 0.6. 

The Mach 0.6 case contains localized regions of supersonic flow as evidenced by tables 
3-4 and by the three-dimensional contour of Mach 1 furnished in figure 2. Nonetheless, the 
statistical properties of the flow remain largely unaffected by compressibility effects. This is 
shown in figures 3-6 which track the time histories of several statistical variables obtained 
from 96 3 DNS. The time histories for skewness (fig. 3), An (fig. 4), and total kinetic energy 
(fig. 5) at Mach numbers 0.0, 0.1, 0.4 and 0.6 are virtually superimposed on each other. 

Flatness and skewness are affected the most by compressibility effects in these simula- 
tions. Figure 3 indicates that the skewness which corresponds to an isotropic turbulent state 
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monotonically increases with Mach number. It is -0.50 at Mo = 0 and has increased to -0.46 
at M 0 = 0.6. Before the flow has reached a state of isotropic turbulence, the time evolu- 
tion of skewness at all Mach numbers are indistinguishable from each other. The physical 
system has equilibrated after approximately one third the total computation time. While 
not reaching an equilibrium value, it is nonetheless worthwhile to point out that the flatness 
parameter decreases by 2% as the Mach number is raised from 0.0 to 0.6 as seen in figure 6. 

Figure 5 illustrates the decay of turbulent kinetic energy (J* cPcc) as a function of 
time. This decay is a natural consequence of viscous damping. After a brief initial increase, 
R\ continuously decreases in time, (fig. 7), with no sign of stabilizing. On the other hand, A 
which is representative of the smaller eddies, increases in time (fig. 4). This indicates that 
energy in the higher wavenumbers is being depleted by the molecular viscosity. 

Tables 3-5 summarize the results of direct simulations of compressible isotropic turbu- 
lence for Mo = 0.1, 0.4 and 0.6, for t=0.1145 on three different grids^. Incompressible results 
are included for comparison. On all the grids, the compressible data converges to the in- 


M 0 

E 

e 

* V max 

Hfr(Sk)) 

(M) 

M max 

0.0 

40.26 

157.4 

0.00 

-0.424 

0.00 

0.00 

0.1 

40.82 

158.2 

0.17 

-0.440 

0.07 

0.21 

0.4 

41.09 

160.4 

1.50 

-0.428 

0.28 

0.84 

0.6 

41.32 

162.3 

3.30 

-0.406 

0.43 

1.26 


Table 3: Summary of direct simulations on a 64 3 grid with initial aver- 
age Mach numbers of 0.0, 0.1, 0.4 and 0.6 at t=0.1145. 


Mo 

E 

e 

* v|max 

li ir ( sk )> 

(M) 

Mmax 

0.0 

40.35 

154.3 

0.00 

-0.506 

0.00 

0.00 

0.1 

40.49 

155.5 

0.14 

-0.505 

0.07 

0.23 

0.4 

40.79 

157.0 

1.17 

-0.493 

0.28 

0.93 

0.6 

41.04 

158.3 

2.82 

-0.477 

0.42 

1.39 


Table 4: Summary of direct simulations on a 96 3 grid with initial aver- 
age Mach numbers of 0.0, 0.1, 0.4 and 0.6 at t=0.1145. 


compressible results as the Mach number is driven towards zero. As expected, the divergence 
of velocity no longer vanishes, and is now an increasing function of M 0 . 

DNS at Mq = 0.6 was not performed on the 128 3 grid. 
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M o 

E 

e 

WBjttBBS 

<i(r(Sk)> 

(M) 


m 

40.29 

153.98 

0.00 

-0.521 

Esa 

0.00 

IB 

40.39 

154.84 

0.14 

-0.518 

ESI 

0.21 


40.78 

156.54 

1.11 

-0.505 

EH 

0.86 


Table 5: Summary of direct simulations on a 128 3 grid with initial 
average Mach numbers of 0.0, 0.1 and 0.4 at t=0.1145 


While the dissipation is approximately the same on the two finer grids, the consistently 
lower values on the coarsest grid confirm the previously stated conclusion that a 64 3 grid 
cannot resolve all the length scales. As a function of increasing Mach number, the trace of 
Sk increases, the dissipation decreases, while the total kinetic energy increases very slightly. 

The results in tables 3-5 are averages over several DNS runs with different initial seeds. 
A given seed uniquely determines the initial velocity distribution, and therefore the pressure 
and temperature fields. Variations of the seed are only necessary to eliminate the statistical 
uncertainty due to the random velocity distribution. The distribution of velocity on two 
different grid sizes are different even when the initial seed is the same. 

Skewness is even more sensitive to the grid refinement than the dissipation as witnessed 
by its decrease from a value of -0.505 to one of -0.521 on 96 3 and 128 3 grids respectively. This 
might be a result of the greater sensitivity of the fluctuating velocity field spatial derivatives 
to slight inaccuracies in the flow variables. 


5.2. Data Analysis 

Using the data generated from the previously discussed DNS of compressible homogeneous 
turbulence at low Mach numbers, the proposed subgrid-scale (SGS) model is now validated. 
Models relate the subgrid-scale stresses - which are not available to a large-eddy simulation 
code - to the large scale velocities which are known. These velocities are simply the Favre- 
filtered velocities introduced earlier. The Favre-filtered velocities are calculated by filtering 
the resolved DNS velocity field with a Gaussian spatial filter of width A = A/Ax/, where 
A Xf is the grid spacing on the fine grid. For convenience, A c and A / refer to the filter width 
A non-dimensionalized with respect to the coarse and fine grid spacing respectively. 

Perturbed velocity fluctuations on the fine grid are the difference between the fully re- 
solved velocity and the filtered ones, given by 

v' = v - v. (70) 

From x 1 and v, subgrid-scale stresses based on DNS, refered to as exact, are calculated. 
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These include the Leonard, cross and Reynolds subgrid-scale stresses given by Eqs. (20)- 
(22). However, these subgrid-scale stresses themselves do not directly affect the evolution of 
the system. The momentum equation is only influenced by the divergence of the subgrid-scale 
stresses (i.e., the vector level). Similarly, v • (V • r) (the scalar level) is a better represen- 
tation of the dissipation terms in the energy equation than are the stresses. Consequently, 
correlations are performed on the tensor, vector and scalar levels. Ideally, high correlations 
are desired on all levels. 

The data analysis proceeds in multiple stages. First, the exact stresses calculated from 
the DNS are injected down to the coarse grid, along with the filtered velocities. The modeled 
subgrid-scale stresses are then calculated (excluding the model constants) on the coarse grid. 
Some variables must be filtered a second time (e.g. cross stress terms). Rather than calculate 
them on the fine grid (which is not available to the large-eddy simulation codes), a Gaussian 
filter is applied to v on the coarse grid with a filter width of A c . Consistency between the 
coarse and fine filter widths is achieved by insuring that 

Nj_ 

A c N c 

where Nf and N c are respectively the number of nodes along one direction of the fine and 
coarse grids. This guarantees that the filtering on the coarse and fine grids is performed 
over the same region in physical space. Derivative evaluations on both the coarse and the 
fine grid are based on Fourier collocation. Calculations by McMillan and Ferziger (1979) 
indicate that the model constants are sensitive to the accuracy of the derivative evaluations. 
A general trend that has been observed is that the Smagorinsky constant is lowered when 
derivative quantities are evaluated more accurately. Our constants are therefore expected to 
lie in the lower range of the values obtained by McMillan (1980). 

Next, the model constant, Cr, is calculated. Unfortunately, the constants can be calcu- 
lated by a wide variety of algorithms, each with its own merits. Moreover, for each algorithm, 
the constants can be evaluated from tensor, vector or scalar information. Therefore, criteria 
must be established to identify the best method. A key test is that L + C should be Galilean 
invariant. To make use of this fact, an additional constant, Cc, is introduced as an extra 
factor in the subgrid-scale cross stress model. A self-consistent method of calculating the 
constants must reproduce Cc = 1 to satisfy the Galilean invariance property stated above 
(Speziale 1985). Additional tests are performed on coarse grids with varying degrees of 
refinement which further decrease the number of choices. A thorough discussion of model 
constants is the subject of the next subsection. 

Once a single or a multiple set of model constants have been determined, the model 
subgrid-scale stresses are calculated and correlated with the exact subgrid-scale stresses 
calculated from the DNS after injection onto the coarse grid. The correlations are performed 
for each type of subgrid-scale stress individually, and for the total stress (L + C + R). Strong 
differences in the correlation coefficients relating total stresses are noticed depending on 
whether or not the Leonard stresses are included. Finally, the correlations obtained from 
the proposed model are compared with the linear combination model, which has been shown 
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to be one of the best models available for incompressible isotropic turbulence. Correlation 
coefficients are calculated based on the two pairs of constants that are obtained from the 
above considerations. The set that is finally retained corresponds to the highest levels of 
correlation of the total stress on the vector and scalar levels. These matters are treated more 
completely in a later subsection. 

To avoid a possible confusion of terminology when referring to variables being compared 
against each other, superscripts m and e are sometimes used. They respectively refer to 
modeled and exact (based on DNS) variables at the tensor, vector and scalar levels. 


5.2.1. Model Constants 

The proposed model given by equations (31) and (42) has two undetermined coefficients. 
The constant, Cr } is associated with the modeled subgrid-scale Reynolds stress, R, while 
Pr T is associated with the thermal heat flux. 

Although the cross stress model has no constant associated with it, it is nonetheless 
multiplied by a constant Cc . This is done in the hope of reducing the number of schemes 
by which the constants can be calculated. A good model should reproduce a cross stress 
constant of one to guarantee Galilean invariance. Once the constants have been determined, 
Cc is set to one and forgotten. Because the flow is isotropic, constants are expected to 
be the same for the three diagonal stress components, the three off-diagonal components 
and the three vector components. Therefore, the values presented in the tables below are 
averaged over the appropriate components. In the tables, D refers to averaged diagonal 
components, OD to averaged off-diagonal components, V to averaged vector components 
and S to averaged scalar components. Similar averagings are performed for the correlation 
coefficients. 

The two constants (Cr and Cc) are calculated using two techniques - each applied on 
the tensor, vector and scalar levels. One method enforces equality of the rms levels of the 
exact and modeled stresses. This is done for each individual subgrid-scale stress (i.e., the 
subgrid-scale Reynolds stress and the cross stress). Hereafter this approach is referred to 
as RMS. The second method utilizes multiple linear least square regression (LSQ) between 
the exact (18) and the modeled (31) total subgrid-scale stresses to determine the constants. 
Table 6 summarizes the results obtained from incompressible data. The three cases presented 
are identical except for the initial random number seed. The constants are independent of 
the detailed velocity statistics. These results are based on a vector level comparison between 
the modeled and exact stresses. Both RMS and LSQ produce Cc close to unity as required. 
Unfortunately this prevents a rational choice from being made between the two approaches. 
A more complete set of LSQ constants is presented in table 7. They are computed at Mach 
numbers of 0.0, 0.1, 0.4 and 0.6 on the tensor, vector and scalar level. Computations are 
performed on a coarse grid of 16 3 . 
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LSQ 

RMS 

Seed 

Cr 

Cc 

Cr 

C c 

i 


1.04 

0.023 

jjfeKl 

2 


1.03 

0.022 


3 


1.03 

0.022 

1.03 | 


Table 6: Model constants calculated by LSQ and RMS between exact 
and modeled total stresses (L + C+R). Results are based on three iden- 
tical incompressible simulations except for the random initial velocity 
distributions. Calculations are on the vector level on a 16 3 grid. 


At first glance, Cc is near unity at both the scalar and vector levels. However whereas 
on the vector level, the constant remains within 4% of unity for all Mach numbers, this is 
not the case on the scalar level where Cc is a decreasing function of Mach number. This 
trend is present in the data generated from both random seeds. Although not presented 
here, the rms cross stress model constant is also near unity when calculated based on vector 
level stresses. Therefore, a preferred method for the determination of the model constants 
is still not possible. 


5.2.2. Filter Width and Grid Coarseness 

Confirmation of the numerical evidence presented by McMillan and Ferziger (1979) that 
A c = 2 is the best filter width is given in table 8 (Mo = 0.1). The criterion used to 
determine the validity of the filter width is that Cc MS must remain close to unity on the 
vector level. Only when A c = 2 is Cc near one. Similar results hold for LSQ constants. The 
constants also vary with respect to the coarse grid on which the LES is to be performed. 
Table 9 summarizes the model LSQ and RMS constants evaluated from modeled stresses 
calculated on 16 3 and 32 3 grids on the vector level. The data ( M 0 = 0.1) shows that Cr 
varies by 30% when the grid size on which the modeled stresses are calculated ranges from 
16 3 to 32 3 . On the other hand, large eddy simulations using finite-difference algorithms 
might be performed on grids as large as 128 3 . Unless a subgrid-scale model is found which 
produces constants independent of the coarse grid size, LES simulations will run the risk 
of producing the wrong results. Perhaps a more complicated dependence of the modeled 
subgrid-scale Reynolds stresses on A is required. 

The best model constants will produce the highest correlations between the modeled and 
exact subgrid-scale stresses on all levels (tensor, vector and scalar). Based on the previous 
discussion, only constants calculated on a vector level are adequate because they produce 
a Cc of unity. Unfortunately, a clear cut choice between RMS and LSQ constants cannot 
be made because Cc (on the vector level) is nearly one in both cases. Rather than make 
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an arbitrary choice, both sets of constants are considered when correlating exact against 
modeled stresses. For reference, the constants used henceforth are 

Cr SQ = 0.012 (72) 

C* MS = 0.023. (73) 



Seed 1 

Seed 2 

Mo 


0.0 

0.1 

0.4 

0.6 

0.0 

0.6 


D 

1.32 

1.32 

1.32 

1.31 

1.32 

1.31 

C c 

V 

1.04 

1.04 

1.02 

1.00 

1.04 

1.00 


s 

1.00 

1.00 

0.95 

0.873 

0.971 

0.934 


D 

0.018 

0.018 

0.018 

0.018 

0.016 

0.015 

c R 

V 

0.012 

0.012 

0.012 

0.012 

0.012 

0.012 


s 

0.015 

0.015 

0.015 

0.014 

0.015 

0.015 


Table 7: LSQ model constants. Filter widths are A / = 12 and A c = 2. 



LSQ 

RMS 

A/ 

6 

12 

f 

6 

12 


A c 

1 

2 

■9 

1 

2 


C R 

0.007 

1 

0.020 

0.019 

0.023 

0.034 

C c 

0.31 

BB 

1.33 

0.82 

1.03 

1.13 


Table 8: LSQ and RMS model coefficients between exact and modeled 
stresses on the vector level at Mo = 0.1. Results are obtained with fine 
filter widths of 6, 12 and 24 while maintaining the proper ratio of 6 
between fine and coarse widths. The coarse grid is 16 3 . 


5.2.3. Correlations 

Correlation coefficients have long been a preferred diagnostic tool for estimating the reliability 
of the modeled stresses. However, the subgrid-scale stresses have been separated into various 
components (L, C, R), each modeled separately and although the correlation of any of these 
components against their models might be excellent, it is still possible for the correlation 
of the exact total stress against the modeled total stress to be less impressive. Such is the 
case, for instance, when two stress components have opposite signs and partially cancel each 
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LSQ 

RMS 

Grid 


EMI 

c R 

C c 

16 3 


1.03 

0.023 

1.03 

32 3 

Biil 

1.03 

0.013 

1.02 


Table 9: LSQ and RMS model constants based on subgrid-scale stresses 
evaluated on 16 3 and 32 3 coarse grids. 


other out. As a final note, before the specific correlation coefficients are presented, one must 
always be attentive to the actual relationship between the exact and the modeled variable, 
even when the correlation coefficient is relatively high (say, above 70%). A correlation 
coefficient as high as 70% may not be as good as it seems. For example, the correlation 
between the functions y == x and y = exp(— x) in the interval [0,1] is approximately 70% 
although they are qualitatively different functions! As a consequence, correlations are only 
deemed good if the correlation coefficient is above, say, the 90% level. 

For convenience, the compressible subgrid-scale model is restated here: 

uRij = 2C R pA 2 (S kl S kl f 2 (Si, - (74) 

iRii = 0 (75) 

Cij = -p(viVj - ViVj ) (76) 


The correlation coefficients presented in table 10 between exact and modeled R and C, 
are independent of the model constants. The correlation coefficients are insensitive to the 
average Mach number variation. The neglect of jRij appears to be a good approximation; the 
direct simulations show it to be several orders of magnitude smaller than the thermodynamic 
pressure. For example, for all of the compressible isotropic turbulence simulations conducted 
in this study, 


(V ■ jR) r 
(Vp)™ 


< 3 x 10 


-3 


(77) 


and, hence, the effect on the isotropic part of the Reynolds stress tensor is dominated by the 
thermodynamic pressure. 


In most of the literature on subgrid-scale models, the Leonard stress has been omitted 
from the total stress correlations on the grounds that it is calculated exactly (Bardina, 
Ferziger, and Reynolds 1983, McMillan 1980, McMillan and Ferziger 1979). However, it has 
recently been shown by Speziale (1985) that the combination L + C is Galilean invariant, 
while the Leonard and cross stresses, individually, are not. Therefore we feel that correlations 
of the total stress should include the Leonard stress. Table 11 summarizes the correlation 
coefficients between the total stress including and excluding the Leonard stress. Results 


20 












Mo 


0.1 

0.4 

0.6 


D 

31 

31 

31 

R 

OD 

26 

26 

25 


V 

22 

22 

22 


S 

45 

45 

45 


D 

89 

89 

89 

C 

OD 

91 

91 

91 


V 

80 

80 

80 


S 

75 

74 

72 


Table 10: Correlations between exact and modeled stresses, R and C, 
at Mach 0.1, 0.4 and 0.6. 


are presented at Mach 0 and Mach 0.6. When the Leonard stresses are left out, correlation 
coefficients similar to those of Bardina are obtained on all levels. However, the inclusion of L 
decreases the correlations at the vector level by approximately 30%. Table 11 substantiates 
that the correlation coefficients (with and without L) are nearly independent of the initial 
average Mach number. 

Correlation coefficients between (C + R) e and various combinations of the modeled 
stresses using constants calculated by LSQ and RMS are summarized in table 12. The 
second column indicates the model against which the total stress is being compared. Al- 
though at first glance RMS based constants perform better at the tensor level when the 
total modeled stress is considered, at the vector and scalar levels, the trend is reversed. 
Correlations at the vector and scalar level are higher by 4% using the LSQ constants. 

When the constants are selected based on LSQ, table 12 indicates that the correlations 
on all levels are highest when all modeled components are included. However, the dynamic 
evolution of the large scale velocities only brings into play the stresses on the vector and 
scalar level. Therefore the coefficients which produce the maximum correlations of total 
stress on these two levels should be chosen. This leads to an optimum choice of 


C R = 0.012 (78) 

calculated by least squares fit of the total stress on the vector level. 

Conversion of Cr to the standard form currently used in incompressible LES ' produces 
a Smagorinsky constant of 0.092. McMillan (1979) obtained a value of Cs = 0.10 when 
spectral collocation derivative computations were employed. This constant corresponds to 
a scalar level evaluation based on RMS. On the vector level, McMillan calculated a higher 

Uln a number of reports (McMillan and Ferziger 1979, Bardina et al. 1983), the subgrid-scale Reynolds 
stress model is proportional to C§. In these cases, the relationship between Cs and Cr is Cr = y/2 Cj. 
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M , o = 

0.0 

M o = 

0.6 | 


L+C+R 

C+R 

L+C+R 

C+R 

D 

93 

82 

93 

81 

OD 

80 

85 

79 

84 

V 

46 

72 

46 

71 

S 

56 

73 

56 

74 


Table 11: Comparison of correlation coefficients of the exact total stress 
versus its model. The modeled terms are computed on a 16 3 coarse grid. 
Each case is presented with and without the inclusion of the Leonard 
subgrid-scale stress terms (calculated exactly). Both the incompressible 
and the M=0.6 case are shown to illustrate the weak influence of Mach 
number on the correlation coefficients with Cr = 0.012. 




Least squares 

RMS 



C R = 0.0122 

Cr = 0.023 

Exact 

Model 

D OD V S 

D OD V S 

C + R 

C 

C + R 

78 82 68 62 

81 84 71 70 

78 82 68 62 

88 83 67 67 


Table 12: Correlation of the exact total stress (C + R) e with various 
models. The modeled stresses are defined in equations (74)-(76). 


Smagorinsky constant of 0.13. This value can be obtained from the present data by using 
x>R instead of £>R However as shown above, the correlations of the total subgrid- 
scale stress would be lower. 

Initial tests of the subgrid-scale heat flux model produced a turbulent Prandtl number 
in the range of 0.4 to 0.5. For this initial study, we take 

Pr T = 0.5 (79) 

which is a value that has been used in previously published large-eddy simulations of turbu- 
lent flows with thermal convection (cf. Eidson 1985). A more accurate calculation of Pr t 
could be accomplished in a compressible flow with significant mean temperature gradients; 
this is beyond the scope of the present study and must await future research. 
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6. Large-Eddy Simulation of Compressible Isotropic Turbulence 


Now, in order to demonstrate the efficacy of the subgrid-scale models derived in this paper, 
a large-eddy simulation of compressible isotropic turbulence is conducted. Since most high- 
speed compressible turbulent flows have significant regions where the turbulence statistics are 
quasi-incompressible, it is important that the models perform well for weakly compressible 
turbulence - the type of flow considered in the last section for the a priori tests. However, 
it is well known that a priori tests only provide a relatively weak gauge for the performance 
of subgrid-scale models in an actual large-eddy simulation (see Hussaini, Speziale and Zang 
1990). It is therefore important to examine their performance in an actual large-eddy sim- 
ulation, particularly for a case where significant compressibility effects are exhibited in the 
turbulence statistics. 

Direct simulations of compressible isotropic turbulence were conducted corresponding to 
the initial conditions Re = 250, (M 2 )\j 2 = 0.1, (/> rmJ )o = 0, (T rms ) o = 0.0626 and two values 
of Xo : 0 and 0.2 which, respectively, correspond to initial values for (.Ra)o of 30.0 and 26.3. 
Here, x = E c /(E 1 + E c ) where E 1 and E c are the incompressible and compressible parts of 
the turbulent kinetic energy, respectively. The direct numerical simulations were performed 
on a 96 3 mesh. Characteristic energy and dissipation spectra associated with the DNS are 
shown in figures 8 (a)-(c) and 9 (a)-(c) for Xo = 0 and xo = 0.2, respectively. Figure 8 (b) 
clearly shows that x remains very small for all time, with only very slight modification of the 
energy spectrum in time. In contrast, figure 9 (b) shows an initial cascade of the spectrum 
towards the smaller scales, followed by a strong energy dissipation at the smaller scales of vf . 
The dissipation spectra ( k 2 E{k') ) shown in figures 8 (a) and 9 (a) demonstrate that both the 
small and large scales are well resolved. Comparing figures 8 (c) and 9(c), the incompressible 
energy E 1 is the same at t = 0 and t = 1.6, but is influenced by compressibility effects at 
the intermediate times. This influence is characterized by a slight decrease in E 1 for Xo > 0. 

Integral properties of these two cases are plotted in figures 10 (a)-(e). These figures 
contrast the two runs corresponding to Xo = 0 and xo = 0.2. The higher compressibility 
has a variety of effects on the flow. The total kinetic energy decays at a slightly slower rate, 
while both the skewness and the flatness are decreased. In other words, finite compressibility 
drives the flow more towards a Gaussian state. Both the integral scale Lj and the Taylor 
microscale An become smaller as Xo increases. Lastly, the decay rate of the microscale 
Reynolds number is slower for finite Xo- A more detailed study of these effects awaits future 
research. It would be particularly useful to compute these statistical quantities based on the 
solenoidal and irrotational components of velocity. 

As noted by Frisch and Orszag (1990), the three dimensional vorticity is formed of 
tubular-like structures. This is graphically represented in figure 11 which shows a volume 
rendering of lj 2 . Volumetric rendering is a visualization technique whereby rays are projected 
from the eye through the flow (which emits and absorbs light). The tubular structures of 
vorticity are contrasted with the spherical structure of (V • v) 2 shown in figure 12. Thi§ 
difference is to be expected since the dilatation satisfies an isotropic wave equation to lead- 
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ing order which shows no preferential direction. On the other hand, the vorticity stretching 
occurs along the axis of the vorticity vector, which then becomes a preferred direction. 

The results for the direct simulation with the strongest compressibility effects (% 0 = 0-2) 
were filtered and injected onto a coarse 32 3 grid for comparison with an LES which was also 
performed on a 32 3 grid using the subgrid-scale models derived herein. In this manner, a 
direct comparison can be made between the results of the LES and the direct simulation 
for a flow where the turbulence statistics exhibit significant compressibility. The following 
turbulence statistics were compared: 


(1) the integrated average of UiUi and 

(2) the integrated average total and isotropic vortex stretching (denoted by UiSijWj and 
u> 2 V • v, respectively), 

(3) the mean turbulence Mach number defined by 

(4) the mean compressible, incompressible and total turbulent kinetic energy denoted by 
E c , E 1 and Etot> respectively, 

(5) the level of compressibility x defined by x = E c /E tot, and 

(6) the rms of the thermodynamic pressure, density, and temperature denoted by P Tmt , 
p rm „ and T rm „ respectively. 

These quantities represent a good choice of turbulence statistics to monitor the effects of 
compressibility. 

In figure 13 (a)-(f), a direct comparison of these statistics for the DNS and LES is made 
for a filter width A/ = 2. It is clear that the LES does an excellent job in reproducing 
the results of the DNS with the possible exception of x ■ It should be noted that x exhibits 
acoustic oscillations and hence is a difficult quantity to predict accurately; nonetheless, the 
LES yields results that are in good qualitative agreement with the direct simulation. The 
most striking result is how well the compressible turbulent kinetic energy and dilatational 
terms are captured. 

It was found that a change in the filter width Af to 1 or 3 - and an adjustment of the 
SGS model constants Cr and Prj of up to 25% - only led to a mild degradation of the 
results. However, the subgrid-scale models did play a crucial role in obtaining the accurate 
results shown in figure 13 (i.e., a 32 3 direct simulation is substantially under-resolved). To 
illustrate this point, the results of a direct simulation on the coarse 32 3 grid (i.e., an LES 
for A / = 0) is shown in figure 14 for the same test case. It is clear from this figure that the 
coarse grid DNS does a poor job in capturing the incompressible as well as the compressible 
turbulence statistics. Consequently, it is the adequate performance of the subgrid-scale 
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models in draining the proper amount of energy from the filtered fields (to account for the 
presence of unresolved scales) which leads to the excellent results obtained in figure 13. 
Considering the degree of compressibility exhibited by the turbulence statistics for the test 
case considered in this section, it would appear that the feasibility of the proposed subgrid- 
scale models has been established. A more extensive parametric study of this subgrid-scale 
model is under investigation by Dahlburg et al. (1990). 


7. Conclusion 


New subgrid-scale models for compressible turbulence have been developed and tested against 
the results of direct numerical simulations of compressible isotropic turbulence. These com- 
pressible subgrid-scale models, which were based on the Favre-filtered equations of motion 
for an ideal gas, contain two dimensionless constants and reduce to the linear combination 
model of Bardina, Ferziger and Reynolds (1983) in the incompressible, isothermal limit. The 
subgrid-scale stress model constant was found to assume the value of Cr — 0.012 which gives 
rise to correlations between the exact and modeled stresses that were above 70% on the ten- 
sor, vector and scalar levels - a correlation which compares favorably with those obtained 
in earlier work on the subgrid-scale modeling of incompressible turbulent flows. Another en- 
couraging feature lies in the fact that these constants and their associated correlations were 
found to be comparatively insensitive to a mean Mach number in the range 0 < M 0 < 0.6. 

The results of a coarse grid 32 3 LES of compressible isotropic turbulence conducted with 
the subgrid-scale models derived in this paper were shown to be in excellent agreement with 
those obtained from a 96 3 direct simulation. These results are extremely encouraging since, 
for the case considered, on average 25% of the turbulent kinetic energy was compressible. 
Furthermore, the ability for the LES to accurately capture the dilatational statistics of the 
flow was quite surprising. 

Future research will be directed on several fronts. The large-eddy simulation of a com- 
pressible turbulent flow with mean temperature gradients could lead to refinements in the 
subgrid-scale heat flux model. Furthermore, the large-eddy simulation of compressible, ho- 
mogeneous shear flow could yield new insights into the performance of these subgrid-scale 
models. Near- wall modifications will also be implemented that allow for the large-eddy sim- 
ulation of compressible, wall- bounded turbulent flows. While further improvements are still 
possible, we believe that the essential foundation for the large-eddy simulation of compress- 
ible turbulent flows has been established in this study. With future research, compressible 
LES could have a profound impact on the analysis of supersonic and hypersonic flows of 
aerodynamic importance. 
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Appendix 


Initial energy spectrum 


Both the incompressible and compressible direct simulations impose a specified energy- 
spectrum on the initial random velocity distribution. For comparative purposes, the data 
was obtained from tabular data found in Comte-Bellot and Corrsin (1971). They tabulate 
the function Eu(k) which is related to the energy spectrum E(k) by 


E(k) = 


i k 3 — 

2 dk 


ldE n 

k dk 


(Al) 


Unfortunately the data is noisy, so a least squares fit is performed on log£ n , expressed as 
a fourth order polynomial in log k. The final form obtained for E X1 is 


l°g(£ n ) = 2.64359 - 0.72602(log k) - 0.32585(log A:) 2 
+0.03525(log A:) 3 - 0.02344(log A:) 4 . 


(A2) 


Calculation of model constants 


Before correlating the total exact subgrid-scale stress against its model, the model con- 
stants must be determined. There are many ways of accomplishing this among which two 
are retained. The total modeled stress is written as a linear combination of modeled terms 
CiT™: 

r m = E CiT™ (A3) 

i= 1 

while the exact total stress is simply 

T ' = ll T l (A4) 

t=l 

The unknown constants to be determined are the C { . The first approach adopted is to 
calculate the root mean square of the pairs C ^ and r? and equate them for each value of 
i. The constants thus take the values 


Ci = 


ov? 

t 

<7 r m 


This method is referred to by RMS in the text. 


(A5) 


A least square method applied to the total stress as a whole is an alternative approach. 
In this case, the norm 

ir-T^HI£W-C,T;)||’ (A6) 

1=1 

is minimized with respect to the coefficients C{. This gives rise to a linear system in the 
coefficients which can be solved by direct methods if the number of constants is not too 
large. For the subgrid-scale model considered in the text, n = 3. 
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Definitions 


For reference purposes, several statistical definitions are provided here. All variables are 
defined on a three-dimensional grid and are subscripted by a single index i for convenience. 
The average of a function Ti is 

< A7 > 

iV i= 1 

As a function of the average, the rms of T is 

r rm , = -{?))’)■ (AS) 

Correlation coefficients are fundamental in evaluating subgrid-scale models. The correlation 
coefficient between two functions .Fand Q is 

C{F,G) = ( A9 ) 

rm«yrmj 


29 


List of Figures 


Fig. 1 Energy spectra at i = 0.0586 for both aliased and de-aliased calculations for 64 3 , 96 3 , 
and 128 3 grids. 

Fig. 2 Three-dimensional Mach 1.0 contours of the DNS at t = 0.0586. The direct simulation 
was performed on a 96 3 grid. 

Fig. 3 Time history of one third the trace of Sk for M 0 of 0.0, 0.1, 0.4 and 0.6. Direct 
simulation was performed on a 96 3 grid. CBC initial conditions. 

Fig. 4 Time history of one third the trace of A^ for M 0 of 0.0, 0.1, 0.4 and 0.6. Direct 
simulation was performed on a 96 3 grid. CBC initial conditions. 

Fig. 5 Time history of kinetic energy for Mq of 0.0, 0.1, 0.4 and 0.6. Direct simulation was 
performed on a 96 3 grid. CBC initial conditions. 

Fig. 6 Time history of one third the trace of FI for M 0 of 0.0, 0.1, 0.4 and 0.6. Direct 
simulation was performed on a 96 3 grid. CBC initial conditions. 

Fig- 7 Time history of R\ for Mq of 0.0, 0.1, 0.4 and 0.6. Direct simulation was performed 
on a 96 3 grid. CBC initial conditions. 

Fig. 8 Energy and dissipation spectra of 96 3 DNS at t = 0, 0.2, 0.4, 0.8, 1 .6. jfc 0 = 10, xo = 0.0, 
(M t ) o = 0.1, ( R x )o = 30.0, Re = 250. 

(a) Dissipation, (b) Compressible energy, (c) Solenoidal energy. 

Fig. 9 Energy and dissipation spectra of 96 3 DNS at t = 0,0.2, 0.4, 0.8, 1.6. Jfc 0 = 10, xo = 0.2, 
(M t ) o .= 0.1, (R\)o = 26.4, Re = 250. 

(a) Dissipation, (b) Compressible energy, (c) Solenoidal energy. 

Fig. 10 Turbulence statistics for 96 3 DNS at t = 0.8. Parameters are k 0 = 10, xo = 0.0 and 
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Fig. 3 Time history of one third the trace of Sk for Mq of 0.0, 0.1, 0.4 and 0.6. Direct 
simulation was performed on a 96 3 grid. CBC initial conditions. 
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Fig. 4 Time history of one third the trace of A ^ for M 0 of 0.0, 0.1, 0.4 and 0.6. Direct 
simulation was performed on a 96 3 grid. CBC initial conditions. 




Fig. 5 Time history of kinetic energy for of 0.0, 0.1, 0.4 and 0.6. Direct simulation was 
performed on a 96 3 grid. CBC initial conditions. 






Fig. 6 Time history of one third the trace of FI for M 0 of 0.0, 0.1, 0.4 and 0.6. Direct 
simulation was performed on a 96 3 grid. CBC initial conditions. 
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Fig. 7 Time history of R\ for Mq of 0.0, 0.1, 0.4 and 0.6. Direct simulation was performed 
on a 96 3 grid. CBC initial conditions. 




Fig. 8 Energy and dissipation spectra of 96 3 DNS at t = 0, 0.2, 0.4, 0.8, 1.6. k 0 = 10, %o = 0.0, 
(Af t ) 0 = 0.1, (i2,\)o = 30.0, Re = 250. 

(a) Dissipation, (b) Compressible energy, (c) Solenoidal energy. 
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Fig. 10 Turbulence statistics for 96 3 DNS at t = 0.8. Parameters are k 0 = 10, xo = 0.0 and 
Xo = 0.2, (Mt) o = 0.1, (prmt)o = 0, (T rmt ) o — 0.0626, Re = 250. 

(a) Etot = Ej C + E 1 , (b) trace of skewness tensor, (c) trace of flatness tensor, (d) trace 
of Taylor length tensor, (e) integral length scale, (f) microscale Reynolds number. 
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Fig. 12 Volumetric representation of (V • v) 2 . Parameters are identical to those of Fig. 10. 
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Fig. 14 Comparison between LES results ( , , ) and results from a 96 3 DNS 

(A, □, O) injected onto the 32 3 LES grid. Filter width A c = Aj = 0. 
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